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Abstract 

Using Mathematica 3.0, the Schrodinger equation for bound states is solved. The method of solution is 
based on a numerical integration procedure together with convexity arguments and the nodal theorem 
for wave functions. The interaction potential has to be spherically symmetric. The solving procedure 
is simply defined as some Mathematica function. The output is the energy eigenvalue and the reduced 
wave function, which is provided as an interpolated function (and can thus be used for the calculation 
of, e.g., moments by using any Mathematica built-in function) as well as plotted automatically. The 
corresponding program schroedinger . nb can be obtained from f ranz . schoeberl@univie . ac . at. 



PACS: 03.65.Ge 
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1 Introduction 

The Schrodinger equation is one of the fundamental equations (of motion) in physics. 
Unfortunately, exact analytic solutions may be found only in exceptional cases. A wide 
area of application has been and still is nonrelativistic potential models, which describe 
bound-state properties of hadrons, considered as bound states of quarks interacting via 
some spherically symmetric potential. Detailed discussions may be found in Refs. [|. 
Fortunately, the two-body Schrodinger equation with spherically symmetric potential 
V(r), where r = |x| and x is the relative coordinate of the constituents, can be reduced 
to an ordinary differential equation for the reduced wave function y n /(r), describing a 
bound state of radial and orbital-angular-momentum quantum numbers n and £, resp. 
In natural units, where h = c = 1, the (nonrelativistic) two-body Schrodinger equation 
in configuration space for the reduced wave function, which is normalized according to 



reads 




; dr [y n/ {r)Y = 1 , (1) 
o 

d 2 £{£ + 1) 



y n A r ) = E n,eynA r ) ■ ( 2 ) 



£ = 0,1,2,... denotes the angular-momentum quantum number, \i the reduced mass, 

_ mi 777.2 

A* = : > 

777i + 7772 

and E n i is the energy eigenvalue, with 77 = 0, 1, . . . counting the number of nodes of the 
bound-state wave function within (0, 00), which corresponds to the radial excitations. 

Here, a Mathematica |§ notebook is constructed for the numerical solution of the 
reduced Schrodinger equation. The solution of differential equations with Mathematica 
consumes more computational time than with Fortran since the Mathematica program 
is not compiled. On the other hand, it is rather complicated to handle graphics within 
Fortran. Moreover, using Fortran one has to write a program for any new potential, to 
compile it, and to link it. Thus, it is much more tedious to study various potentials with 
Fortran than with Mathematica. In Mathematica, it is rather easy to define a function, 
to calculate matrix elements, and to use graphics tools. 

The outline of this paper is to discuss first the method of finding the solutions, then 
to demonstrate the application of Mathematica, and, finally, to give, in the Appendix, 
the program listing. (It may be obtained as Mathematica notebook schroedinger .nb 
from f ranz . schoeberl@univie . ac . at.) The method applied here is based on Ref. [3J]. 

2 Method of Solution 

Rewriting Eq. (@) in more convenient form, the differential equation to be integrated is 

VnA r ) = \V&( r ) - £ nA Vn,i(r) , (3) 

with the effective potential 

V eS (r) = 2^V(r) + £ -^^- 
and the scaled energy eigenvalue 

£ n ,e = 2/7 E n> e . 
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For the energy eigenvalue E to be bounded from below, the potential V(r) has to be, 
in any case, less singular than — 1/r 2 . Moreover, for potentials V(r) such that r 2 V(r) is 
analytic (which implies that V(r) is less singular than 1/r 2 ), the general (normalizable) 
solution of the differential equation (|3|) is given as a power-series expansion of the form 

y(r) oc r m [1 + 0(r)] . 

The asymptotic behaviour of the not normalized reduced wave function y{r) for r — ► 
is therefore 

lim y(r) = r i+1 . 

The normalization condition (|]) forces the reduced wave function y(r) to approach 
zero when r — > oo. Consequently, the main idea of finding the energy eigenvalues E and 
corresponding reduced wave functions y is to perform a systematic scan for increasing 
values of e in Eq. ([3]), looking for those values of e which allow y{r) — > for r — > oo. 
Without loss of generality, near the origin y may be assumed to be positive: y(0+) > 0. 

• For e low enough, V e s (r) — e is certainly positive and thus y(r) — ► +oo for r — > oo. 

• Increasing e, the divergence of y(r) for r — > oo will weaken. 

• Increasing e, V e g(r)— £ will become negative in certain regions of r. If this region is 
large enough, it may happen that y(r) vanishes for r — ► oo (cf. Fig. [I]): the lowest 
bound-state energy E 0> g for given £ has been found. 
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Figure 1: Qualitative behaviour of potential V e s(r) and reduced ground-state wave function y(r) 



• Increasing e, y(r) will cross somewhere and behave like y(r) — > — oo for r — ► oo. 

• Increasing e further, it may happen that y(r) approaches from below for r — ► oo: 
the energy eigenvalue Eij of the first radial excitation for given £ has been found. 

In more detail our procedure works as follows: In principle, the integration of Eq. (|^) 
starts at the origin. For singular potentials, however, it has to start at some value r = S 
close to but different from the origin and to respect, of course, the boundary conditions 

y(S) = 5 m , 



= (£+1)5' 



(4) 
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First, note that there is a classical turning point such that V e g(r) > e n ^ for all r > r c \, 
the classical turning point r c \ is the largest value of r solving the equation V e s (r) = e n ,i- 
Secondly, note that 

sign y"(r) = sign y(r) for all r > r c i . 

This means, the reduced wave function 2/ is convex (concave) if it is positive (negative). 
Thus, at some point r> > r c i, y(r>) > and y'(r>) > implies y(r) — > +00 for r — > 00 
while y(r>) < and £/(r>) < implies j/(r) — > —00 for r — > 00. Clearly, for both cases 
the integration can be stopped. 

The level of excitation found in the course of the integration procedure is identified, 
according to the well-known nodal theorem, by the number n of nodes of the (reduced) 
radial wave function: the ground-state wave function has no node at all (n = 0), while 
the wave function of some radial excitation has a finite number of nodes (n = 1, 2, . . .). 

In order to locate a desired bound state, define an interval for the energy e, by fixing 
appropriate lower and upper bounds E L and E^, resp. The routine starts to integrate, 
with the Runge-Kutta method || , the differential equation (H) at r = 5, respecting the 
boundary conditions (|4]) and using, for the energy e, the arithmetic mean (E L + Eu)/2. 
It counts the number of nodes, n, detected within a prescribed interval (E^, E\j) during 
the integration and changes E^ and E\j appropriately in an iterative procedure: if these 
bounds are chosen badly, e converges to that bound which lies next to the true energy; 
in order to cover nevertheless the desired bound state, the corresponding bound has to 
be changed]] 



3 Applications 

The notebook is called schroedinger . nb. It identifies r c \ by determination of the local 
minimum of the effective potential at largest r. For numerical reasons, it uses an upper 
bound on this minimum by adding one stepsize h to it. Of course, for some potentials, 
like all pure power-law potentials, this local minimum may be determined analytically. 
For instance, for some potential of the form V(r) = r k , k G N, the minimum resides at 

1(1+1) ] V(fc+2) 
k fj, 

For more complicated potentials, the local minimum has to be determined numerically. 
This is done by the module xwmill; the arguments of the latter are orbital excitation £, 
stepsize h of the numerical integration, as well as stepsize we it and starting value xrat 
of the minimum search. Both I and h are taken over from the calling procedure schroe, 
which calls xwmill with stepsize weit=0.5 and guessed minimum at xrat=20. (These 
values are rather good starting values for the minimum search; they may be changed in 
the module schroe.) 

Upon starting the notebook schroedinger . nb, one is provided with explanations of 
the usage as well as a description of the definitions used. The procedure goes as follows: 

Input: schroedinger .nb. 

Output: description of usage and definitions. 

Input: potential vl [r_] : = . . . ; e.g., for the harmonic-oscillator potential vl [r_] : =r~2. 



1 The program listing in the Appendix gives more details. 



4 



Input: schroe [0 , 20 , 2 , 1 , . 01 , 1 , 1] . The arguments of this routine are: 

lower bound on the energy eigenvalue e in chosen units of energy (= 0), 

upper bound E\j on the energy eigenvalue e in chosen units of energy (= 20), 

number n of nodes = radial excitations (= 2), 

angular momentum i = orbital excitations (=1), 

stepsize h in chosen units of energy (= 0.01), 

masses m 1; m 2 in chosen units of energy (mi = m 2 = 1). 

Output: energy eigenvalue E n> i in chosen units of energy, etc. 

The above example generates as output the seeked numerical results and — if desired— 
a plot of the (not normalized) reduced wave function y (here called yschr) as function 
of r (here called x): 

E = 12.99999714, L = 1, N = 2, 

Integration steps = 589, h = 0.01, del =0.001, el = 0, eu = 20, 

Largest x, upper integration limit, XMAX = 5.881, 

Smallest x, lower integration limit, XMIN = 0.001. 

The not normalized reduced wave function is yschr [x] . 

The normalization factor is given by: 

1/NIntegrate [yschr [x] "2 , {x,del ,xmax}] 




We are now able to use Mathematica built-in functions to evaluate matrix elements 
like (1/r). 

Input: NIntegrate [1/r yschr [r] "2 ,{r ,del ,xmax}] / 
NIntegrate [yschr [r] "2 , {r , del , xmax}] . 

Output: 0.625982. 

Table p] lists results of a test run for the harmonic-oscillator potential while Table |2] 
shows results for the Coulomb potential for varying stepsizes h and starting points del. 

Obviously, y(0) has to vanish for all states while y'(0) is nonvanishing only for I = 
states — which provides a (trivial) additional consistency check. Of course, the required 
computational time depends strongly on the desired accuracy. 

In summary, the Mathematica notebook developed here provides an easy-to-handle 
procedure for computing energies and wave functions of (nonrelativistic) bound states. 
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Table 1: Energy eigenvalues, moments, as well as wave function and its derivative at the origin for the 
harmonic-oscillator potential V(r) = r 2 . The values chosen for the parameters are m\=mi = \ GeV, 
h = 0.01 GeV -1 , E L = GeV, and E v = 20 GeV. The exact energies are E n<e = 4n+2£+3 GeV; the 
discrepancies of the computed values only reflect the /i-dependent accuracy of our numerical method. 
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2.999997 
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1.9977 


0.0000 


2.2567 
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7.000003 


0.9403 


1.9966 


0.0000 


3.3851 
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10.999999 


0.8369 


1.9958 


0.0000 


4.2314 
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4.999995 


0.7523 


0.6667 


0.0000 


0.0000 


1 


1 


9.000001 


0.6770 


0.6667 


0.0000 


0.0000 


2 


1 


12.999997 


0.6260 


0.6667 


0.0000 


0.0000 





2 


7.000003 


0.6018 


0.4000 


0.0000 


0.0000 


1 


2 


10.999999 


0.5588 


0.4000 


0.0000 


0.0000 


2 


2 


15.000005 


0.5266 


0.4000 


0.0000 


0.0000 



Table 2: Energy eigenvalues, moments, and derivative of the wave function at the origin for the IS, 2S, 
and 2P states of the Coulomb potential V(r) = —1/r as function of the stepsize h. The starting point 
6 is calculated from 5 = h/10; the numerical values chosen for the parameters are mi = = 1 GeV, 
£l = —1 GeV, and E\j = +1 GeV. The corresponding exact results are given in parentheses. 



State 


IS 




2S 




2P 




h [GeV -1 ] 


0.1 


0.05 


0.1 


0.05 


0.1 


0.05 


-E n ,t [GeV] 


0.249973 


0.249996 


0.062496 


0.062496 


0.062504 


0.062504 


(0.25) 




(0.0625) 




(0.0625) 




(1/r) [GeV] 


0.4999 


0.4999 


0.1250 


0.1250 


0.1250 


0.1250 


(0.5) 




(0.125) 




(0.125) 




(1/r 2 ) [GeV 2 ] 


0.4947 


0.4974 


0.06185 


0.06221 


0.02082 


0.02082 


(0.5) 




(0.0625) 




(0.02083) 




[y'(0)] 2 [GeV 3 ] 


0.4897 


0.4949 


0.06123 


0.06190 


2.10~ 6 


4.10- 7 


(0.5) 




(0.0625) 




(0) 
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A Program Listing 

(* 

Author: Franz F. Schoberl 

Institute for Theoretical Physics 
University of Vienna 
Boltzmanng. 5 
A- 1090 Vienna 

E-MAIL : Franz . Schoeberl@Univie . ac . at 

SOLVING THE SCHRODINGER EQUATION FOR BOUND 
STATES WITH MATHEMATICA 3.0 

This software is based on: 

" Solving the Schrodingcr Equation for Bound States" by 
P. Falkcnstciner, H. Grosse, Franz F. Schoberl, and P. Hertel 
Computer Physics Communication 34 (1985) 287-293. 

The energy and the reduced radial wave function for a bound state with given number of nodes nO 
and angular momentum 1 is calculated. The potential has to be spherically symmetric. 

Usage: 

el lower bound of the energy, 

eu upper bound of the energy, 

nO radial excitations, number of nodes, 

1 orbital excitations, 

feh ....error on the energy, built in as feh=0. 00001, 

h integration stepsize, determines the number of the integration steps 

and thus the accuracy of the determined energy 

ml, m2 constituent masses, 

schroe[el,eu,n0,l,h,ml,m2] calling procedure 

The output reduced wave function (not normalized) is called yschr[x]. 
The plot of the wave function is called yschrplot. 
Example, harmonic oscillator: 

(1) Define the potential: vl[x_]:= x~ 2 

(2) Start solving the Schrodinger equation: schroe[-l, 20, 1,2, 0.1, 0.336, 0.336] 

The above procedure will solve the equation for the first radial and second orbital excitation. 
The equation to be solved is (with h = c=l, fx = ™*™^ 2 ): 

+ ^7^) + V{ r )]ynA r )i yschr[x] is the not normalized y n ,e(r), 

^n,e,m(r) = *^Y t , m {il), / (y n A.r)) 2 dr = 1. 

Jo 

The accuracy can be increased by decreasing h, this increases the number of integration steps. 
The higher the number of integration steps the more accurate the eigenvalues as well as the eigen 
functions. The reduced wave function will be plotted in addition. A measure for the accuracy is also 
the shape of the wave function. It should vanish for the largest numerical x, otherwise one has to 
decrease h. 

If you run schroe[el,eu,n0,l,h,ml,m2] you automatically will be asked if you like to plot the reduced 
wave function. 

*) 
(* 

xwmill calculates the minimum of the potential most to the right. The minimum 

is called xwmin. xrat is some guessed x- value most to the right 

(here I have used x-rat = 20). wcil is the stepsize of the minimum search 

(here I have used weil = 0.5). 

*) 

xwmill [11_, hl_, wcil_, xratl_, wwl_] := 

Module[{l = 11, h = hi, weit = weil, xrat = xratl, ww = wwl}, 



del = h/10; 
xs = xrat; 

If[xs < del, Goto [11]]; 
If[l < 1, Goto [3]]; 
Label[lj; 

If[xs<weit + del, xs = del; Goto [2]]; 

ms = ww*vl[xs] + 1*(1 + l)/xs' 2; 

rs = ww*vl[xs + weit] + 1*(1 + l)/(xs + weit)' 2; 

Is = ww*vl[xs - weit] + 1*(1 + l)/(xs - weit)" 2; 

If[rs>ms && ms>ls, xs = xs - weit; Goto[l]]; 

If[rs<ms && ls>ms, xs = xs + weit; Goto[l]]; 

If[rs = Is , Goto [2]]; 

If[rs>ms && ls>ms, If[wcit<h, Goto [2]]; 
weit = weit/10; Goto[l]]; 

Print ["SOMETHING IS WRONG, ANALYZE POTENTIAL"]; 
Goto[10]; 
Label [3]; 
Label [1 a]; 

If[xs<wcit + del, xs = del; Goto [2]]; 

ms = ww*vl[xs]; 

rs = ww*vl[xs + weit]; 

Is = ww*vl[xs - weit]; 

If[rs>ms && ms>ls, xs = xs - weit; Goto[l a]]; 
If[rs<ms && ls>ms, xs = xs + weit; Goto[l a]]; 
If[rs = Is , Goto [2]]; 

If[rs>ms && ls>ms, If[wcit<h, Goto [2]]; 
weit = weit/10; Goto[l a]]; 

Print ["SOMETHING IS WRONG, ANALYZE POTENTIAL"]; 

Goto[10]; 

Label[ll]; 

Print ["XMIN TOO SMALL"]; 

Goto[10]; 

Label [2]; 

xwmin = xs; 

Label[10]; 

Return["end"]], 

(* 

The program schroe 
*) 

schroe[ell_, eul_, n01_, 11_, hl_, ml., m2_]:= 

Module[{el = ell*ww, eu = eul*ww, nO = nOl, 1 = 11, h = hi, wl 
w2 = m2}, 

(* 

Determining the minimum of the potential most to the right. 

*) 

ww = 2*wl*w2/(wl + w2); 
xwmill[l, h, 0.5, 20, ww]; 
xwmin = xwmin + h; 
del = h/10; 
fch = 0.00001*ww; 

(* 

Defining dim = Vcff - c, y" = (Vcff - e) y. 

*) 

diffl[xxx_, L, ccn_] := ww*vl[xxx] + 1*(1 + l)/xxx~ 2 - een; 
Label [300]; 
sch = eu - el; 
cps = (el + eu)/2; 



Starting the integration with the boundaries y(del) = del" (1 + 1), 
y'(del) = (1 + l)der 1, and with the energy eps equal to the arithmetic 
mean of the preceding lower and upper, bounds of energy el and eu. 

x = del; 

y = x'(l+ 1); 

yp =(1 + l)*x~l; 

yold = 1; 

nOx = 0; 

If the desired accuracy (prescribed error feh) has been obtained, the bound 
state energy is taken as the arithmetic mean of the last el and eu. 

If[seh<fch, Goto[l]]; 

Integrating y" = (Vcff - e)y one step h further with the Rungc - Kutta 
method 

Label [2]; 

al = yp*h; bl = difR[x, 1, cps]*h*y; 

a2 = (yp + bl/2)*h; hh = diffl[x + h/2, 1, eps]*h; 

b2 = hh*(y + al/2); a3 = (yp + b2/2)*h; b3 = hh*(y + a2/2); 

a4 = (yp + b3)*h; 

x = x + h; u2 = diffl[x, 1, eps]; b4 = u2*h*(y + a3); 
y = y + (al + 2*a2 + 2*a3 + a4)/6; 
yp = yp + (bl + 2*b2 + 2*b3 + b4)/6; 

Counting the number of nodes by nOx until the prescribed nO is reached. 

If[y*yold>0, Goto[3]]; 
nOx = nOx + 1; 
If[n0x>n0, Goto[4]]; 
Label [3]; 
yold = y; 

If the following condition is not fullfilled, x is greater then the classical 
turning point (u2 has the value of Veff - eps at the point x) . 

If[(u2<0 || x<xwmin), Goto [2]]; 

If (after stating that x is greater than the classical turning point) y*yp is 
greater than (i.e., y and yp have the same sign), one is sure that y goes to 
infinity, without having additional nodes. Otherwise one has to integrate 
further, 

z = y*yp; 

lf[z<0, Goto[2]]; 

If y goes to infinity, a new el is established by eps. 

el = eps; 
Goto [300]; 

If nox exceeds nO, a new eu is established by eps. 



Label [4]; 
eu = eps; 



Goto [300]; 



In the following lines the wave function y is calculated using the above 
calculated bound state energy (the last eps) by the same method as above. In 
addition y is stored in feldl at x which is stored in xcoord, and the number 
of integration steps is counted by j . 

Label[l]; 
cp = eps; 
j = 0; 
Label [20]; 

j = j + 1; 
Mdl [j] - y; 

xcoord[j] = del + (j - l)*h; 

ji = j; 

xmax = xcoord[jl]; 

Integrating y" = (Veff - e)y one step h further with the 
Runge - Kutta method 

al = yp*h; bl = diffl[x, 1, cps]*h*y; 

a2 =(yp + bl/2)*h; hh = diffl[x + h/2, 1, eps]*h; 

b2 = hh*(y + al/2); a3 = (yp + b2/2)*h; b3 = hh*(y + a2/2); 

a4 =(yp + b3)*h; 

x = x + h; u2 = diffl[x, 1, eps]; b4 = u2*h*(y + a3); 
y = y + (al + 2*a2 + 2*a3 + a4)/6; 
yp = y + (bl + 2*b2 + 2*b3 + b4)/6; 
lf[y*yold>0, Goto[30]]; 
nOx = nOx + 1; 
If[n0x>n0, Goto[40]]; 
Label [30]; 
yold = y; 

If[(u2<0 || x<xwmin), Goto[20]]; 

z = y*yp; 

lf[z<0, Goto[20]]; 
Label [40]; 

The reduced radial wave function yschr obtained from the interpolation 
of the data is stored in feldl and in xcoord 

yschr = Interpolation[Table[{xcoord[j],feldl[j]}, {j, 1, jl}]]; 
xmax = xcoord [jl]; 

Output of the resulting energy eigenvalue and the input data. 
Print [ 

StylcForm["E = ", 

FontColor -> RGBColor[0.996109, 0, 0], Font Weight -> "Bold", 
FontSize -> 16], 
StyleForm[N[ep/ww, 10], 

FontColor -> RGBColor[0.996109, 0, 0], Font Weight -> "Bold", 
FontSize -> 16], 
StyleForm["L = ", 

FontColor -> RGBColor[0.996109, 0, 0], Font Weight -> "Bold", 

FontSize -> 16], 

StyleForm[l, 

FontColor -> RGBColor[0.996109, 0, 0], Font Weight -> "Bold", 
FontSize -> 16], 



StyleForm["N = ", 

FontColor -> RGBColor[0.996109, 0, 0], Font Weight -> "Bold", 

FontSize -> 16], 

StylcForm[nO, 

FontColor -> RGBColor[0.996109, 0, 0], FontWeight -> "Bold", 
FontSize -> 16], 

Integrationstcps = ", jl, "h = ", h, " del =", 

del, " el = ", ell, " eu = ", eul, 

StyleForm[" Largest x, upper integration limit, XMAX = ", 

FontColor -> RGBColor[0, 0, 0.996109], FontWeight -> "Bold", 

FontSize -> 16], 

S ty lcForm [xmax , 

FontColor -> RGBColor[0, 0, 0.996109], FontWeight -> "Bold", 
FontSize -> 16], 

StylcForm[" Smallest x, lower integration limit, XMIN = del = ", 
FontColor -> RGBColor[0, 0, 0.996109], FontWeight -> "Bold", 
FontSize -> 16], 
StyleForm[del, 

FontColor -> RGBColor[0, 0, 0.996109], FontWeight -> "Bold", 
FontSize -> 16],".", 

StyleForm["The reduced not normalized wave function is yschr[x]. 

The normalizationfactor is given by: 

1 /NIntegrate[yschr[x] * 2,{x,del,xmax}]" , 

FontColor -> RGBColor[0.996109, 0.500008, 0.250004], 

FontWeight -> "Bold", FontSize -> 16]]; 

Preparing the plot of the reduced wave function yschr[x] 

zz = InputString["You like to plot the (not normalized) reduced wave 
function? Type yes and click OK, otherwise click just OK "]; 
If[zz == "yes", yschrplot — Plot[yschr[x], {x, del, xmax}, 
PlotStyle -> GrayLevel[0], AxesLabel -> {"x", "yschr"}, 
DcfaultFont -> {"Times-Bold", 12}, 
Background -> RGBColor[0.996109, 0.996109, 0], 
PlotLabel -> FontForm["Not normalized 
reduced wave function", {"Helvetica-Bold", 14}]], 

Print [StyleForm[" OK NO PLOT", FontWeight -> "Bold", FontSize -> 18, 
FontColor -> RGBColor[0, 0.500008, 0]]]]; 
Return []] 



